{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variational Quantum Linear Solver\n",
    "*Copyright (c) 2022 Institute for Quantum Computing, Baidu Inc. All Rights Reserved.*\n",
    "## Background\n",
    "\n",
    "System of linear equations is a basic yet extremely useful tool in mathematics. An example is that in economics, you can model the economy using linear equations. Also, it provides simple estimation for large system of non-linear systems. Hence solving system of linear equations is an important task.\n",
    "\n",
    "Variational Quantum Linear Solver (VQLS) is a variational quantum algorithm for solving system of linear equations. It's a classical-quantum hybrid algorithm that can run on recent Noisy Intermediate-Scale Quantum (NISQ) devices. To be more specific, given a matrix $A$ and a vector $\\boldsymbol{b}$, our goal is to find a vector $\\boldsymbol{x}$ so that $A \\boldsymbol{x} = \\boldsymbol{b}$. Using VQLS, we can obtain a quantum state that is proportional to $\\boldsymbol{x}$, i.e. a normalised vector $|x\\rangle = \\frac{\\boldsymbol{x}}{\\lVert \\boldsymbol{x} \\rVert_2}$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Principle\n",
    "\n",
    "Solving linear equations in the quantum setting is different to the general setting due to the requirement in quantum computing that we can only apply unitary operators to a quantum state. For the input matrix $A$, we need to decompose it to a linear combination of unitary operators $A = \\sum_n c_n A_n$ where each $A_n$ is a unitary operator. For the input vector $\\boldsymbol{b}$, we need to assume that it's a quantum state that can be prepared by unitary operator $U$, i.e. $U|0\\rangle = |b\\rangle$.\n",
    "\n",
    "![VQLS](vqls.png)\n",
    "\n",
    "We can see that the algorithm consists of two parts. On a quantum computer, we prepare a parameterized quantum circuit (PQC) $V(\\alpha)$ and compute the loss function $C(\\alpha)$, then on a classical computer, we minimize parameters $\\alpha$ until the loss function is below a certain threshold, denoted as $\\gamma$. At the end, we output the target quantum state $|x\\rangle$. The main idea behind the algorithm is that PQC $V(\\alpha)$ gives us a quantum state $|\\psi(\\alpha)\\rangle$, circuit $F(A)$ then computes how similar $A|\\psi(\\alpha)\\rangle$ and $|b\\rangle$ are, which is what the loss function $C(\\alpha)$ is measuring. When the loss is small, $A|\\psi(\\alpha)\\rangle$ and $|b\\rangle$ are very close, it means $|\\psi(\\alpha)\\rangle$ and the target $|x\\rangle$ are very close, so we output $|\\psi(\\alpha)\\rangle$ as an approximation to $|x\\rangle$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Paddle Quantum Implementation\n",
    "\n",
    "We use the `Circuit` class in Paddle Quantum and optimizer in Paddle Paddle to implement VQLS. For the quantum part, we use built-in `complex_entangled_layer` ansatz to build our PQC $V(\\alpha)$. To compute the loss function, we use Hadamard Test and Hadamard-Overlap Test which utilizes `oracle` gate to implement the controlled-$A_n$ gates. For the classical optimization part, we used Adam optimizer to minimize the loss function.\n",
    "\n",
    "User can use toml file to specify the input to the algorithm, matrix $A$ and vector $\\boldsymbol{b}$, stored as '.npy' files. You can run the following code to randomly generate a $n\\times n$ matrix $A$ and vector $\\boldsymbol{b}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Here is a randomly generated A:\n",
      "[[4.1702199e+00+7.203245j   1.1437482e-03+3.0233257j\n",
      "  1.4675589e+00+0.9233859j  1.8626021e+00+3.4556072j\n",
      "  3.9676747e+00+5.3881674j ]\n",
      " [2.0445225e+00+8.781175j   2.7387592e-01+6.704675j\n",
      "  4.1730480e+00+5.5868983j  1.4038694e+00+1.9810148j\n",
      "  8.0074453e+00+9.682616j  ]\n",
      " [8.7638912e+00+8.946067j   8.5044211e-01+0.39054784j\n",
      "  1.6983042e+00+8.781425j   9.8346835e-01+4.2110763j\n",
      "  9.5788956e+00+5.3316526j ]\n",
      " [6.8650093e+00+8.346256j   1.8288277e-01+7.5014434j\n",
      "  9.8886108e+00+7.4816566j  2.8044400e+00+7.892793j\n",
      "  1.0322601e+00+4.4789352j ]\n",
      " [2.8777535e+00+1.3002857j  1.9366957e-01+6.7883554j\n",
      "  2.1162813e+00+2.6554666j  4.9157314e+00+0.5336254j\n",
      "  5.7411761e+00+1.4672858j ]]\n",
      "Here is a randomly generated b:\n",
      "[4.191945 +6.852195j  3.1342418+6.9232264j 6.9187713+3.1551564j\n",
      " 9.085955 +2.9361415j 5.8930554+6.9975834j]\n"
     ]
    }
   ],
   "source": [
    "n = 5\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "np.random.seed(1)\n",
    "A = np.zeros([n, n], dtype=\"complex64\")\n",
    "b = np.zeros(n, dtype=\"complex64\")\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        x = np.random.rand() * 10\n",
    "        y = np.random.rand() * 10\n",
    "        A[i][j] = complex(x, y)\n",
    "    x = np.random.rand() * 10\n",
    "    y = np.random.rand() * 10\n",
    "    b[i] = complex(x, y)\n",
    "np.save(\"./A.npy\", A)\n",
    "np.save(\"./b.npy\", b)\n",
    "print(\"Here is a randomly generated A:\")\n",
    "print(A)\n",
    "print(\"Here is a randomly generated b:\")\n",
    "print(b)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "User can specify the parameters of the VQLS in the toml file. They are `depth`, `iterations`, `LR` and `gamma`, which correspond to the number of layer in the PQC $V(\\alpha)$, number of iterations of the optimizer, learning rate of the optimizer and threshold of the loss function to end optimization early. By entering `python vqls.py --config config.toml` one could solve the linear system. Here we present an example of an online demo. First, define the content of the configuration file as follows, user can try out different settings by changing the parameters of `test_toml`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_toml = r\"\"\"\n",
    "# The path of the input matrix A. It should be a .npy file.\n",
    "A_dir = './A.npy'\n",
    "# The path of the input vector b. It should be a .npy file.\n",
    "b_dir = './b.npy'\n",
    "# The depth of the quantum ansatz circuit.\n",
    "depth = 4\n",
    "# Number optimization cycles.\n",
    "iterations = 200\n",
    "# The learning rate of the optimizer.\n",
    "LR = 0.1\n",
    "# Threshold for loss value to end optimization early, default is 0.\n",
    "gamma = 0\n",
    "\"\"\""
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we run the VQLS:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 88%|████████▊ | 176/200 [02:03<00:16,  1.43it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Threshold value gamma reached, ending optimization\n",
      "Here is x that solves Ax=b: [ 1.3475237 -0.7860472j   0.22970617-0.88826376j -0.35111237-0.31225887j\n",
      "  0.07606918+1.2138402j  -0.729564  +0.48393282j]\n",
      "This is actual b: [4.191945 +6.852195j  3.1342418+6.9232264j 6.9187713+3.1551564j\n",
      " 9.085955 +2.9361415j 5.8930554+6.9975834j]\n",
      "This is Ax using estimated x: [4.185339 +6.8523855j 3.1297188+6.923625j  6.924285 +3.1467872j\n",
      " 9.092921 +2.932943j  5.8879805+6.999589j ]\n",
      "Relative error:  0.0008446976\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "import argparse\n",
    "import os\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "os.environ[\"PROTOCOL_BUFFERS_PYTHON_IMPLEMENTATION\"] = \"python\"\n",
    "\n",
    "import toml\n",
    "import numpy as np\n",
    "import paddle\n",
    "from paddle_quantum.data_analysis.vqls import compute\n",
    "\n",
    "paddle.seed(0)\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    config = toml.loads(test_toml)\n",
    "    A_dir = config.pop(\"A_dir\")\n",
    "    A = np.load(A_dir)\n",
    "    b_dir = config.pop(\"b_dir\")\n",
    "    b = np.load(b_dir)\n",
    "    result = compute(A, b, **config)\n",
    "\n",
    "    print(\"Here is x that solves Ax=b:\", result)\n",
    "    print(\"This is actual b:\", b)\n",
    "    print(\"This is Ax using estimated x:\", np.matmul(A, result))\n",
    "    relative_error = np.linalg.norm(b - np.matmul(A, result)) / np.linalg.norm(b)\n",
    "    print(\"Relative error: \", relative_error)\n",
    "    np.save(\"./answer.npy\", result)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Citation\n",
    "\n",
    "```\n",
    "@misc{bravo-prieto2020variational,\n",
    "  title = {Variational {{Quantum Linear Solver}}},\n",
    "  author = {{Bravo-Prieto}, Carlos and LaRose, Ryan and Cerezo, M. and Subasi, Yigit and Cincio, Lukasz and Coles, Patrick J.},\n",
    "  year = {2020},\n",
    "  month = jun,\n",
    "  number = {arXiv:1909.05820},\n",
    "  eprint = {1909.05820},\n",
    "  eprinttype = {arxiv},\n",
    "  doi = {10.48550/arXiv.1909.05820}\n",
    "}\n",
    "```\n",
    "\n",
    "## References\n",
    "\n",
    "[1] “Variational Quantum Linear Solver: A Hybrid Algorithm for Linear Systems.” Carlos Bravo-Prieto, Ryan LaRose, Marco Cerezo, Yigit Subasi, Lukasz Cincio, Patrick J. Coles. arXiv:1909.05820, 2019."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pq-dev",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.15 (default, Nov 10 2022, 13:17:42) \n[Clang 14.0.6 ]"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "5fea01cac43c34394d065c23bb8c1e536fdb97a765a18633fd0c4eb359001810"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
